home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / dassl / ddaini.f next >
Text File  |  1996-07-19  |  7KB  |  258 lines
  1.       SUBROUTINE DDAINI (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR,
  2.      +   IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)
  3. C***BEGIN PROLOGUE  DDAINI
  4. C***SUBSIDIARY
  5. C***PURPOSE  Initialization routine for DDASSL.
  6. C***LIBRARY   SLATEC (DASSL)
  7. C***TYPE      DOUBLE PRECISION (SDAINI-S, DDAINI-D)
  8. C***AUTHOR  PETZOLD, LINDA R., (LLNL)
  9. C***DESCRIPTION
  10. C-----------------------------------------------------------------
  11. C     DDAINI TAKES ONE STEP OF SIZE H OR SMALLER
  12. C     WITH THE BACKWARD EULER METHOD, TO
  13. C     FIND YPRIME.  X AND Y ARE UPDATED TO BE CONSISTENT WITH THE
  14. C     NEW STEP.  A MODIFIED DAMPED NEWTON ITERATION IS USED TO
  15. C     SOLVE THE CORRECTOR ITERATION.
  16. C
  17. C     THE INITIAL GUESS FOR YPRIME IS USED IN THE
  18. C     PREDICTION, AND IN FORMING THE ITERATION
  19. C     MATRIX, BUT IS NOT INVOLVED IN THE
  20. C     ERROR TEST. THIS MAY HAVE TROUBLE
  21. C     CONVERGING IF THE INITIAL GUESS IS NO
  22. C     GOOD, OR IF G(X,Y,YPRIME) DEPENDS
  23. C     NONLINEARLY ON YPRIME.
  24. C
  25. C     THE PARAMETERS REPRESENT:
  26. C     X --         INDEPENDENT VARIABLE
  27. C     Y --         SOLUTION VECTOR AT X
  28. C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
  29. C     NEQ --       NUMBER OF EQUATIONS
  30. C     H --         STEPSIZE. IMDER MAY USE A STEPSIZE
  31. C                  SMALLER THAN H.
  32. C     WT --        VECTOR OF WEIGHTS FOR ERROR
  33. C                  CRITERION
  34. C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS
  35. C                  IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
  36. C                  IDID=-12 -- DDAINI FAILED TO FIND YPRIME
  37. C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS
  38. C                  THAT ARE NOT ALTERED BY DDAINI
  39. C     PHI --       WORK SPACE FOR DDAINI
  40. C     DELTA,E --   WORK SPACE FOR DDAINI
  41. C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
  42. C                  MATRIX INFORMATION
  43. C
  44. C-----------------------------------------------------------------
  45. C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV
  46. C***REVISION HISTORY  (YYMMDD)
  47. C   830315  DATE WRITTEN
  48. C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
  49. C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
  50. C   901026  Added explicit declarations for all variables and minor
  51. C           cosmetic changes to prologue.  (FNF)
  52. C   901030  Minor corrections to declarations.  (FNF)
  53. C***END PROLOGUE  DDAINI
  54. C
  55.       INTEGER  NEQ, IDID, IPAR(*), IWM(*), NONNEG, NTEMP
  56.       DOUBLE PRECISION
  57.      *   X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
  58.      *   E(*), WM(*), HMIN, UROUND
  59.       EXTERNAL  RES, JAC
  60. C
  61.       EXTERNAL  DDAJAC, DDANRM, DDASLV
  62.       DOUBLE PRECISION  DDANRM
  63. C
  64.       INTEGER  I, IER, IRES, JCALC, LNJE, LNRE, M, MAXIT, MJAC, NCF,
  65.      *   NEF, NSF
  66.       DOUBLE PRECISION
  67.      *   CJ, DAMP, DELNRM, ERR, OLDNRM, R, RATE, S, XOLD, YNORM
  68.       LOGICAL  CONVGD
  69. C
  70.       PARAMETER (LNRE=12)
  71.       PARAMETER (LNJE=13)
  72. C
  73.       DATA MAXIT/10/,MJAC/5/
  74.       DATA DAMP/0.75D0/
  75. C
  76. C
  77. C---------------------------------------------------
  78. C     BLOCK 1.
  79. C     INITIALIZATIONS.
  80. C---------------------------------------------------
  81. C
  82. C***FIRST EXECUTABLE STATEMENT  DDAINI
  83.       IDID=1
  84.       NEF=0
  85.       NCF=0
  86.       NSF=0
  87.       XOLD=X
  88.       YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR)
  89. C
  90. C     SAVE Y AND YPRIME IN PHI
  91.       DO 100 I=1,NEQ
  92.          PHI(I,1)=Y(I)
  93. 100      PHI(I,2)=YPRIME(I)
  94. C
  95. C
  96. C----------------------------------------------------
  97. C     BLOCK 2.
  98. C     DO ONE BACKWARD EULER STEP.
  99. C----------------------------------------------------
  100. C
  101. C     SET UP FOR START OF CORRECTOR ITERATION
  102. 200   CJ=1.0D0/H
  103.       X=X+H
  104. C
  105. C     PREDICT SOLUTION AND DERIVATIVE
  106.       DO 250 I=1,NEQ
  107. 250     Y(I)=Y(I)+H*YPRIME(I)
  108. C
  109.       JCALC=-1
  110.       M=0
  111.       CONVGD=.TRUE.
  112. C
  113. C
  114. C     CORRECTOR LOOP.
  115. 300   IWM(LNRE)=IWM(LNRE)+1
  116.       IRES=0
  117. C
  118.       CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
  119.       IF (IRES.LT.0) GO TO 430
  120. C
  121. C
  122. C     EVALUATE THE ITERATION MATRIX
  123.       IF (JCALC.NE.-1) GO TO 310
  124.       IWM(LNJE)=IWM(LNJE)+1
  125.       JCALC=0
  126.       CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
  127.      *   IER,WT,E,WM,IWM,RES,IRES,
  128.      *   UROUND,JAC,RPAR,IPAR,NTEMP)
  129. C
  130.       S=1000000.D0
  131.       IF (IRES.LT.0) GO TO 430
  132.       IF (IER.NE.0) GO TO 430
  133.       NSF=0
  134. C
  135. C
  136. C
  137. C     MULTIPLY RESIDUAL BY DAMPING FACTOR
  138. 310   CONTINUE
  139.       DO 320 I=1,NEQ
  140. 320      DELTA(I)=DELTA(I)*DAMP
  141. C
  142. C     COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
  143. C     STORE THE CORRECTION IN DELTA
  144. C
  145.       CALL DDASLV(NEQ,DELTA,WM,IWM)
  146. C
  147. C     UPDATE Y AND YPRIME
  148.       DO 330 I=1,NEQ
  149.          Y(I)=Y(I)-DELTA(I)
  150. 330      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
  151. C
  152. C     TEST FOR CONVERGENCE OF THE ITERATION.
  153. C
  154.       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
  155.       IF (DELNRM.LE.100.D0*UROUND*YNORM)
  156.      *   GO TO 400
  157. C
  158.       IF (M.GT.0) GO TO 340
  159.          OLDNRM=DELNRM
  160.          GO TO 350
  161. C
  162. 340   RATE=(DELNRM/OLDNRM)**(1.0D0/M)
  163.       IF (RATE.GT.0.90D0) GO TO 430
  164.       S=RATE/(1.0D0-RATE)
  165. C
  166. 350   IF (S*DELNRM .LE. 0.33D0) GO TO 400
  167. C
  168. C
  169. C     THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
  170. C     M AND AND TEST WHETHER THE MAXIMUM
  171. C     NUMBER OF ITERATIONS HAVE BEEN TRIED.
  172. C     EVERY MJAC ITERATIONS, GET A NEW
  173. C     ITERATION MATRIX.
  174. C
  175.       M=M+1
  176.       IF (M.GE.MAXIT) GO TO 430
  177. C
  178.       IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1
  179.       GO TO 300
  180. C
  181. C
  182. C     THE ITERATION HAS CONVERGED.
  183. C     CHECK NONNEGATIVITY CONSTRAINTS
  184. 400   IF (NONNEG.EQ.0) GO TO 450
  185.       DO 410 I=1,NEQ
  186. 410      DELTA(I)=MIN(Y(I),0.0D0)
  187. C
  188.       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
  189.       IF (DELNRM.GT.0.33D0) GO TO 430
  190. C
  191.       DO 420 I=1,NEQ
  192.          Y(I)=Y(I)-DELTA(I)
  193. 420      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
  194.       GO TO 450
  195. C
  196. C
  197. C     EXITS FROM CORRECTOR LOOP.
  198. 430   CONVGD=.FALSE.
  199. 450   IF (.NOT.CONVGD) GO TO 600
  200. C
  201. C
  202. C
  203. C-----------------------------------------------------
  204. C     BLOCK 3.
  205. C     THE CORRECTOR ITERATION CONVERGED.
  206. C     DO ERROR TEST.
  207. C-----------------------------------------------------
  208. C
  209.       DO 510 I=1,NEQ
  210. 510      E(I)=Y(I)-PHI(I,1)
  211.       ERR=DDANRM(NEQ,E,WT,RPAR,IPAR)
  212. C
  213.       IF (ERR.LE.1.0D0) RETURN
  214. C
  215. C
  216. C
  217. C--------------------------------------------------------
  218. C     BLOCK 4.
  219. C     THE BACKWARD EULER STEP FAILED. RESTORE X, Y
  220. C     AND YPRIME TO THEIR ORIGINAL VALUES.
  221. C     REDUCE STEPSIZE AND TRY AGAIN, IF
  222. C     POSSIBLE.
  223. C---------------------------------------------------------
  224. C
  225. 600   CONTINUE
  226.       X = XOLD
  227.       DO 610 I=1,NEQ
  228.          Y(I)=PHI(I,1)
  229. 610      YPRIME(I)=PHI(I,2)
  230. C
  231.       IF (CONVGD) GO TO 640
  232.       IF (IER.EQ.0) GO TO 620
  233.          NSF=NSF+1
  234.          H=H*0.25D0
  235.          IF (NSF.LT.3.AND.ABS(H).GE.HMIN) GO TO 690
  236.          IDID=-12
  237.          RETURN
  238. 620   IF (IRES.GT.-2) GO TO 630
  239.          IDID=-12
  240.          RETURN
  241. 630   NCF=NCF+1
  242.       H=H*0.25D0
  243.       IF (NCF.LT.10.AND.ABS(H).GE.HMIN) GO TO 690
  244.          IDID=-12
  245.          RETURN
  246. C
  247. 640   NEF=NEF+1
  248.       R=0.90D0/(2.0D0*ERR+0.0001D0)
  249.       R=MAX(0.1D0,MIN(0.5D0,R))
  250.       H=H*R
  251.       IF (ABS(H).GE.HMIN.AND.NEF.LT.10) GO TO 690
  252.          IDID=-12
  253.          RETURN
  254. 690      GO TO 200
  255. C
  256. C-------------END OF SUBROUTINE DDAINI----------------------
  257.       END
  258.